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A  REACTIVE  MULTI-MATERIAL  FIRST  ORDER  NUMERICAL  MODEL  FOR 
HEAT  CONDUCTION  IN  CYLINDRICAL  COORDINATES 


J.A.  Waschl 


ABSTRACT 


A  reactive  multi-material,  first  order  numerical  model  (Known  as 
HEAT)  for  heat  conduction  in  cylindrical  coordinates  has  been  developed  and 
tested.  The  model  allows  the  inclusion  of  a  maximum  of  three  separate 
materials . 


A  second  order  Taylor  series  expansion  was  employed  In  the 
development  of  the  finite  difference  scheme  within  each  single  material 
layer.  Matching  conditons  on  the  interfaces  are  determined  to  first  order. 

The  model  employs  a  zero  order  Arrhenius  burn  function  to  describe 
the  heat  generated  by  the  tnermal  decomposition  of  the  reactive  components. 

It  also  takes  into  account  the  latent  heat  of  fusion  for  each  material. 

External  heat  is  introduced  via  a  constant  temperature  boundary 
condition  or  a  time  dependent  boundary  condition.  The  model  will  permit 
restarts,  although  only  in  the  cases  were  total  decomposition  of  the  reactive 
fraction  has  not  occurred. 
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A  REACTIVE  MULT  I -MATER IAL  FIRST  ORDER  NUMERICAL  MODEL  FOR 


HEAT  CONDPCTIOH  IK  CYLINDRICAL  COORDINATES 


1.  INTRODUCTIOH 


There  are  a  number  of  occasions  during  either  storage  or  service 
life,  when  the  transient  heating  of  a  munition  may  occur  for  any  one  of  a 
number  of  reasons.  These  reasons  include  an  external  fire,  thermal  conduction 
of  heat  from  a  hot  gun  barrel  into  a  loaded  round,  or  even  frictional  heat 
generated  as  air  flows  past  a  missile  moving  at  supersonic  speeds. 

For  this  report  the  heating  of  a  shell,  trapped  within  a  not  gun 
barrel  was  investigated.  As  the  heat  flows  from  the  environment  into  the 
shell,  the  temperature  rises  and  the  confined  explosive  will  decompose  more 
rapidly  producing  further  heat,  if  the  heat  produced  by  this  chemical 
decomposition  cannot  be  transferred  to  the  surroundings  as  rapidly  as  it  is 
produced,  thereby  stabilising  the  temperature  at  an  elevated  value,  the  mass 
of  explosive  will  self-lgnlte. 


The  above  may  occur  if  the  gun  barrel  temperature  is  very  high.  The 
temperature  of  the  explosive  above  which  a  thermal  explosion  is  produced,  is 
known  as  the  critical  temperature.  Given  the  chemical  and  physical  parameters 
of  a  particular  system  the  critical  temperature  may  be  exceeded. 


A  reactive  multi-material,  first  order  numerical  model  (called 
HEAT),  has  been  developed  to  solve  the  two  dimensional  heat  conduction  problem 
for  a  cylindrical  geometry.  This  model  employs  a  second  order  Taylor  series 
expansion  of  the  transient  heat  eguatlon  within  each  single  material  layer, 
continuity  across  material  interfaces  is  achieved  using  a  first  order 
approximation.  HEAT  may  be  employed  to  determine  whether  a  particular  event 
should  be  deemed  critical  or  under  control. 

In  what  follows,  the  model  used  for  the  calculations  is 
formulated.  This  is  followed  by  a  description  of  how  the  HEAT  code  was 
validated.  Further  improvements,  the  results  obtained  and  a  discussion  of 
them  are  then  given. 


2. 


THE  MODEL 


HEAT  was  developed  to  study  the  transient  heating  from  a  hot  gun 
barrel  into  a  failed  to  fire  shell  although  many  problems  of  similar  geometry 
may  also  be  studied.  It  employs  the  cylindrical  geometry  of  a  shell,  but 
neglects  the  heat  sink  action  of  any  attached  fuze  or  warhead.  The 
longitudinal  axis  of  the  shell  is  taken  as  the  axis  of  symmetry,  or  z-axls, 
while  the  radial  direction  is  denoted  as  the  r-axls.  The  metal  and  liner 
thicknesses  (materials  1  and  2,  respectively,  of  Fig.  l)  need  not  be  uniform, 
although  the  high  explosive  must  have  a  constant  radius.  These  assumptions 
were  considered  acceptable  for  the  problem  at  hand. 


in  Fig  l,  an  illustration  of  the  model  in  its  most  general  form  is 
provided,  i.e.  a  cross  section  through  the  cylinder  with  the  axis  of  symmetry 
being  the  left  hand  boundary.  A  maximum  of  three  distinct  material  layers  has 
been  incorporated  into  this  model,  although  this  can  easily  be  increased  by 
modifying  the  code. 


The  partial  differential  equation  describing  unsteady  thermal 
conduction  with  an  internal  heat  source  is  as  follows: 


p 


C  , 

p  au 
i  at 


v2u 


s 


(1) 


where, 

»  -  density  of  the  material 

u  -  absolute  temperature 
cp  -  specific  heat 
lH  -  thermal  conductivity 


and  where  s  is  the  source  term  described  by  the  Arrhenius  burn  function  111: 
S  -  pO  (l  -  «)n  Ze'E/RU  ,  for  r,z  within  explosive 


-  o 


otherwise 


where  n  -  the  order  of  reaction 

<  -  the  fraction  of  explosive  decomposed. 
Q  -  heat  of  reaction 
Z  -  collision  number 
E  -  activation  energy 
R  -  gas  constant 
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If  we  neglect  the  source  term  for  now,  Eqn.  (l)  can  be  written  in 
cylindrical  coordinates  as: 


r  2  2  2 

p  au  3u  i_  au  i__  au  au 
i  at  -  ar2  +  7  ar  +  r2  a(2  +  az2- 


with  cylindrical  symmetry  of  the  solution,  — -  -  o  and  Eqn.  (3) 
simplifies  to:  as 

C  2  2 

p  au  a  u  l  au  a  u 

*r  57  -  71 +  7  57 +  rr  -  r  * 0  • 

3r  a? 


Eqn.  (4)  is  valid  at  all  interior  modes  of  the  model.  At  r  -  o,  the 
physical  solution  is  the  one  that  is  non-singular.  Using  a  power  series 
expansion  for  the  radial  solution  to  equation  (4)  we  get. 


u(r)  -  u  +  u.r  +  u_  r 
ol  2 


For  small  r,  i.e.  4r 


u( 4r )  -  u  +  u, Ir  +  u„»r  + 
o  l  2 


and  a  Taylor  series  expansion  near  zero  gives, 


,  .  2  ,2 
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ar 


Equating  nice  terms  of  Eqns  (5a)  and  (5b): 


u  -  u( o) 

°  3u 
U1  -  57  'o 


i  a2u  . 

U  m  -  -  I 

2  2^2  0 
dr 


Now  Eqn.  (4)  may  be  rewritten  as# 
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Our  power  series  solution  for  r  must  satisfy  this  equation.  Substituting  the 
power  series  for  u  and  obtaining  the  derivative  with  respect  to  r  gives, 


i  f-  (r  _  4u  +.  o(r!  as  r  -»  o 
r  3r  3r  2 
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and  Eqn.  (6)  is  recast  as: 


C  2  2 

p  3u  „  3  u  3  U 
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for  heat  flow  at  the  axis  of  symmetry. 

A  solution  was  found  to  Eqns.  (4)  and  (7)  subject  to  the  boundary 
conditions: 


u(R,z,t)  -  Tj 

u( r, 0, t )  -  Tj,  t>0 

u(r,L,t)  -  Tj 


together  with  the  initial  condition, 


u( r , z, o )  -  T. 


oSr<R,  0<z<L. 


These  conditions  are  illustrated  in  Fig.  1. 


The  numerical  solution  was  obtained  by  first  writing  Eqns  (4)  and 
(7)  in  a  second  order  finite  difference  form  12),  based  on  a  Taylor 


I 


expansion.  The  Alternating  Direction  Implicit  (ADI)  method  13,4,5,6,7,8,9) 
was  then  employed  to  solve  these  finite  difference  equations. 


The  ADI  method  was  employed  due  to  its  unconditional  stability 
(3,5,7,91  and  because  it  is  a  very  time  efficient  method  for  the  solution  of  a 
multidimensional  problem.  ADI  effectively  reformulates  the  finite  difference 
equations  into  an  algebraic  problem  consisting  of  a  set  of  linear  equations 
whose  coefficients  form  a  trldlagonal  matrix.  The  method  involves  solving 
this  set  of  equations  in  each  coordinate  direction  in  turn,  such  that  two  net 
sweeps  occur  in  a  two  dimensional  problem  such  as  this.  The  first  set  would 
be  written  with  one  direction  expressed  implicitly  and  the  other  expressed 
explicitly.  In  the  second  set  the  situation  would  be  reversed.  The  results 
after  these  two  sweeps  is  a  single  time  step  advance. 


The  rz  plane  of  the  model  was  divided  up  into  a  rectangular  grid  of 
Increment  Jr  and  Jz,  in  r  and  z  directions  respectively*.  The  material 
properties  were  defined  at  each  grid  point  or  node  in  the  plane.  Where 
material  Interfaces  occurred,  control  volumes  (5,6)  were  placed  around  the 
boundary  node  of  Interest  (see  Fig.  2).  For  the  control  volume,  the  sum  of 
the  heat  flow  into  the  node  was  equated  to  the  mass  of  the  control  volume 
times  the  specific  heat  capacity  times  the  rate  of  temperature  change 
simulating  a  non-explosive  material.  This  is  first  order  accurate 
differencing.  Second  order  differencing  at  the  interface  may  be  employed  in  a 
later  version  of  HEAT.  For  the  node  shown  in  Fig.  2  this  would  lead  to  the 
following  set  of  first  order  difference  equations  for  r  expressed  lmpllcity 
and  z  expressed  explicitly  and  without  the  source  term. 


untl/J  (-i  A  ) 
1-13  Jr  1 


n+l/2  S  S_  _  2V 

13  Jr  1  Jr  2  It 


♦  un+1/2  <-£  A  ) 
1+13  Jr  2 


n  / 

Ui 3-1  Ti  V 


n  ,  *2  ,  X1  .  2V . 
U13  Jz  Xl  *  Jz  A2  '  it* 


ui3+i  tz  V 


where 


Aj  -  2«  (  i  -  3/2)  Jr Jz 
A2  -  2»  (1  -  1/2)  JrlZ 

V  -  »lr2 JZ( 1  -  5/4)  +  «Jr2JZ  (i  -  3/4) 
1^  -  thermal  coefficient  m  material  1 
u?  -  temperature  at  node  1,3  at  time  n. 


The  values  for  r  and  z  were  determined  after  a  series  of  test  3obs.  it 
was  then  decided  that  a  minimum  of  10  cells  per  material  layer  was 
required. 


z 


Axis  of  symmetry 


Figure  2.  Description  of  the  control  volume  at  <l,j). 


with  the  boundary  conditions  and  initial  conditions  set,  Eqns.  <  sa ) 
and  (8b),  the  tridiagonal  coefficient  matrix  was  solved  employing  the  Thomas 
Algorithm  1 10,11 1.  This  is  a  fast  solution  technique,  in  which  the  structure 
of  a  tridiagonal  matrix  is  used  to  advantage. 


At  the  end  of  each  incremental  time  step,  the  temperature  increase 
due  to  the  source  term  was  determined  at  all  the  nodes  within  the  explosive. 
Given  this  new  temperature  distribution,  melting  was  determined  and,  if 
necessary,  accounted  for,  by  using  the  latent  heat  of  fusion  in  conjunction 
with  the  melting  point  for  each  material.  If,  on  the  other  hand,  the 
temperature  within  the  explosive  had  risen  to  a  very  high  value  (arbitrarily 
set  at  773°K> ,  total  decomposition  was  assumed  and  the  simulation  stopped. 
The  final  time  step  information  was  then  saved. 


It  should  be  noted  that  although  this  code  was  written  as  second 
order  accurate  within  a  single  material,  the  Interface  calculations  are  only 
first  order  accurate.  The  final  model  is  therefore  only  accurate  to  first 
order . 


3 .  VALIDATION 


HEAT  was  validated  by  using  three  separate  techniques.  Each  one  of 
these  techniques  verified  a  particular  aspect  of  the  code.  The  combined 
results  of  these  tests  provided  a  high  level  of  confidence  in  the  model 
predictions . 

Fig.  3  shows  the  geometric  details  of  a  simplified  shell,  in  all 
cases  studied  the  time  step  was  initially  set  to  0 . l  seconds.  If  after  50 
cycles,  the  temperature  rise  was  less  than  10%  at  every  node,  the  time  was 
increased  by  10%.  The  initial  time  step  was  chosen  as  a  compromise  between 
keeping  the  truncation  errors  low  and  still  maintaining  computational  speed. 
This  procedure  allows  the  time  step  to  increase  as  long  as  the  solution  is 
slowly  varying.  The  shell  was  discretlsed  into  a  square  grid  of  0.002  m  (2 
mm)  for  »r  and  iz  giving  each  material  a  minimum  of  lo  cells.  All  cells  must 
be  the  same  size  in  this  model. 


The  first  method  of  validation  Involved  employing  the  model  to  solve 
the  heat  conduction  equation  where  the  solution  was  either  known  or  easily 
recognizable.  This  is  a  standard  technique.  For  Instance,  this  model 
predicted  a  constant  temperature  when  the  initial  and  boundary  conditions  were 
set  equal  at  time  t  -  o  and  the  heat  source  was  inactive.  Conversely,  if  the 
heat  source  was  active,  the  temperature  within  the  explosive  region  of  the 
model  gradually  Increased. 


Problems  which  had  a  particular  symmetry  in  a  plane  perpendicular  to 
the  rz  plane,  displayed  symmetry  in  that  plane  on  output,  as  expected. 
Incorporating  lower  thermal  conductivity  constants  slowed  the  rate  of 
temperature  rise  within  the  model.  All  these  qualitative  observations  provide 
confidence  in  the  model. 


s 


Figure  3.  Simplified  version  of  a  filled  shell  with  a  driving  band. 
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TABLE  1 


Temperatures  at  the  given  radii  and  times  within  a  steel  cylinder  as 
calculated  by  an  analytic  one  dimensional  solution  and  from  a  HEAT 
simulation  in  brackets.  HEAT  valves  are  taken  from  a  cross-sectional  slice  at 
the  centre  line  i.e.  z  =  0.127  m.  Both  results  are  for  a  boundary  temperature 
of  798K  and  an  internal  temperature  of  298K.  The  steel  cylinder  had  a  radius 
of  0.066  m  and  a  height  of  0.254  m. 


RADIUS ( m) 


Time  (s) 

0 

0 . 02 

0.038 

0 . 056 

23.2 

302.3 

321 . 3 

410.6 

630.5 

(  299) 

(309.2) 

(  383.4  ) 

(601.3) 

57.2 

366 . 8 

409.2 

520.1 

689.7 

(313.5) 

(339. 8) 

(  448.2) 

(647.7) 

106  .  9 

481  . 1 

516.0 

601.5 

722.9 

(344.5) 

(372.5) 

(480.4) 

(664.2) 

179.7 

600.5 

622.5 

675.9 

751 . 4 

(386.6) 

(412.2) 

(510.5) 

(677.1) 

286 . 4 

699.4 

710.4 

737.1 

774.8 

(435.6) 

(458.2) 

(544.8) 

(691.5) 
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TABLE  2 


Table  showing  the  time  to  detonation  as  determined  by  an  effectively  l- 
dimensional  calculation  (14)  and  HEAT  using  stated  constant  boundary 
conditions.  The  l-dimenslonal  code  used  5  cells  (0.012  m  each)  In  the  radial 
direction  and  50  cells  (variable  dimensions  between  0.037  m  and  0.002  m)  In 
longitudinal  direction,  and  with  a  time  step  of  0.1  s.  Our  model  used  a 
constant  cell  dimension  of  0.002  m,  33  cells  along  the  radial  direction  and 
127  cells  In  the  longitudinal  direction.  The  time  step  was  monotonlcally 
Increased  from  0.1  s. 


Boundary  Temperature 
(K) 

Time  to  Detonation 
(S) 

1-dlmensional  calculation  2-dlmenslonal  calculation 

(HEAT) 

573 

42.1 

598 

33.3 

35.4 

608 

31.3 

623 

28.9 

32.1 

648 

25.8 

29.1 

698 

24.7 

798 

16.9 

19.2 

2000 

7.0 

8.3 

In  the  second  method,  the  temperature  calculated  from  an  analytical 
solution  for  a  one  dimensional  inert  single  material  was  compared  with  this 
model  in  an  inert  single  material  phase  only.  The  results  are  given  in 
Table  1 . 


A  definite  difference  between  the  analytic  and  numeric  solutions  can 
be  seen  in  Table  1.  This  discrepancy  is  due  mainly  to  the  heat  sinking  effect 
of  the  neighbouring  material  cells  in  the  two  dimensional  case.  Initially  the 
difference  is  not  great  (about  5*),  as  the  heat  conduction  is  predominantly  in 
the  radial  direction.  Later  on,  longitudinal  heat  flow  becomes  important, 
and  the  discrepancy  between  the  one  dimensional  and  the  two  dimensional  models 
grows  to  around  60%.  At  this  stage  the  effect  of  a  two  dimensional 
calculation  is  apparent.  The  result  is  a  lower  calculated  temperature  at  the 
nodes  for  the  two  dimensional  model.  This  demorstrates  the  need  for  a  two 
dimensional  analysis  when  considering  the  temperature  rise  within  a  cylinder. 


In  the  third  method,  HEAT  is  compared  with  the  solution  from  another 
numerical  model  (14!.  Although  the  latter  model  is  two  dimensional,  it 
contains  the  wrong  geometry  for  our  problem  and  hence  can  only  be  U3ed  to 
provide  a  one  dimensional  solution.  Nevertheless,  the  comparison  is  still 
instructive.  Table  2  shows  the  times  to  detonation  for  HEAT  and  for  this  one¬ 
dimensional  model. 


TABLE  3 


Table  showing  the  time  to  detonation  as  determined  by  HEAT.  Either  a  constant 
boundary  condition  or  a  slowly  decreasing  temperature  profile  for  the  boundary 
condition  were  used.  Appendix  A  provides  the  form  of  this  temperature 
profile . 


Initial  Boundary  Temperature 
<K> 

Time  to 

Constant  conditions 

Detonation  (s) 

Decreasing  Conditions 

573 

42.1 

598 

35 . 4 

- 

€  0  8 

57.2 

£23 

32.1 

45.7 

648 

29.1 

35.4 

£98 

24.7 

27.6 

798 

19.2 

20.5 

2000 

8.3 

16.6 

The  time  to  detonation  is  larger  when  calculated  by  HEAT.  This  is 
not  surprising  considering  the  interpretations  made  from  Table  1.  The  effect 
of  including  the  second  dimension  is  to  slow  the  heating  rate  and  hence 
increase  the  time  to  detonation,  as  the  detonation  time  is  the  critical 
parameter  required  from  these  calculations,  the  superiority  of  the  two 
dimensional  model  becomes  apparent.  Even  though  this  time  difference,  as 
shown  here,  is  small,  the  effect  is  real  and  should  be  included  in  the 
calculation. 
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4.  EMBELLISHMENTS .  RESULTS  AMD  DISCUSSION 


A  constant  boundary  temperature  is  unrealistic.  Heat  transfer  from 
the  barrel  to  the  shell,  for  instance,  is  most  likely  to  occur  through  the 
intimate  contact  of  the  driving  band.  Hence  the  boundary  condition  should 
specifically  Include  this  region,  and  for  completeness,  should  also  contain  a 
cooling  function  for  the  temperature  profile  since  the  barrel  will  cool  as  the 
shell  heats  (see  Appendix  A). 

The  cooling  function  was  determined  from  experimental  data  and  gives 
temperature  as  a  function  of  time.  The  length  of  the  driving  band  was 
determined  from  drawings.  The  effect  of  these  changes  can  be  seen  in  Table 
3.  Here,  a  comparison  between  constant  temperature  and  reducing  temperature 
boundary  conditions  is  made.  It  should  be  noted  that  with  the  constant 
boundary  temperature  condition  a  detonation  will  always  occur  provided  that 
this  initial  boundary  temperature  is  greater  than  the  self-ignition 
temperature.  In  the  reducing  temperature  condition,  detonation  only  occurs 
above  a  critical  initial  temperature,  below  this  temperature  self-sustained 
decomposition  cannot  begin.  The  critical  ignition  temperature  in  this  case  is 
a  function  of  the  geometry  of  the  shell  and  the  charge  and  their  thermal 
properties,  as  well  as  of  the  decomposition  properties  of  the  explosive,  Eqn. 
(2). 

In  addition,  the  form  of  the  cooling  function,  the  length  of  the 
driving  band,  and  the  geometry  of  the  problem,  are  all  critical  parameters  in 
determining  whether  a  mass  detonation  will  occur.  For  this  reason,  the 
results  of  Table  3  can  only  be  taken  as  a  guide.  Nevertheless,  the  indication 
is  that  a  significant  change  in  the  predicted  stability  of  the  shell  is  found 
if  the  boundary  conditions  are  more  realistic. 

The  length  of  the  driving  band  in  the  previous  calculations  was 
0.06  m.  As  the  actual  heat  input  may  well  occur  beyond  the  boundaries  of  the 
driving  band,  it  was  decided  to  investigate  time  to  detonation  for  a  longer 
driving  band,  we  chose  an  extended  length  of  o.os  m.  Some  simulations  from 
Table  3  were  then  repeated  using  the  Increased  length.  The  results  are 
compared  in  Table  4 . 


TABLE  4 


Time  to  detonation  as  determined  by  HEAT  for  a  driving  band  of  length 
or  o.oa  m.  The  Initial  boundary  temperatures  are  given.  The  rate  of 
temperature  decrease  is  the  same  for  both  calculations. 


0.06m 

boundary 


Initial  Boundary  Temperature 
(K> 

Time  to  Detonation 
(s) 

Driving  Band  0.06  m  Driving 

Band  o.oa  m 

5  58 

- 

_ 

598 

- 

47.4 

608 

57.2 

€23 

45 . 7 

37.0 

€98 

27.6 

26.1 

798 

20.5 

19.2 

with  the  larger  driving  band,  more  energy  can  be  put  into  the 
system.  As  the  thermal  characteristics  of  the  shell  remain  constant,  the  time 
to  detonation  must  decrease.  It  is  also  interesting  to  note  that  the 
difference  in  the  predicted  times  is  not  large  when  the  initial  temperature  is 
far  above  the  cook-off  temperature.  This  happens  because  at  those 
temperatures  cook-off  occurs  so  rapidly  that  not  all  the  available  energy  is 
used  or  needed.  At  lower  temperatures,  energy  balance  becomes  critical  more 
quickly  for  the  narrower  band  and  hence  the  time  differences  Increase.  The 
extent  of  influence  of  the  driving  band  is  significant  in  the  determination  of 
the  time  to  detonation  and  should  therefore  be  carefully  considered. 


To  improve  efficiency,  the  input  and  output  of  HEAT  was  modified.  A 
useful  feature  that  was  included  allows  the  user  to  restart  the  calculation 
from  the  last  time  step  saved.  This  means  that  the  model  can  be  run  for  a 
small  number  of  cycles  and  the  output  viewed.  The  simulation  can  then  be 
continued  again  if  necessary,  therefore  eliminating  the  need  to  attempt  one 
large  simulation  run  only,  which  is  not  always  desirable. 


5 .  CONCLUSION 


A  model  to  solve  the  two-dimensional  transient  heat  equation  in 
cylindrical  coordinates  has  been  developed  and  programmed  in  FORTRAN  77  to  run 
on  a  VAX  11/760  computer.  This  model  (heat)  has  been  tested  and  compared  with 
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results  for  single  material  or  mul tlmaterlal  media  with  distinct  boundaries. 
The  model  has  been  demonstrated  to  solve  problems  with  either  constant 
boundary  conditions,  or  with  time  dependent  temperature  profiles  along  a 
region  of  the  boundary. 


The  model  uses  second  order  finite  differences  within  a  single 
material  but  only  first  order  finite  differences  at  material  boundaries. 
These  approximations  will  Influence  mesh  size  determination.  They  do  not 
affect  accuracy  of  the  results  directly. 

The  model  should  prove  to  be  valuable  In  the  determination  of  heat 
conduction  Into  filled  shells,  particularly  where  heat  sinks  are  of  concern 
and  where  simpler  one  dimensional  codes  will  not  suffice.  In  these  cases, 
the  calculated  temperature  rises  will  be  more  realistic.  The  model  can  also 
be  used  as  a  basis  for  other  more  general  problems  of  cylindrical  symmetry. 
For  example,  voids  and  tapered  profiles  may  be  incorporated. 
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